% include the BC folder
s = cd; s = [s, '/BC']; path(s, path);
s = cd; s = [s, '/EFG']; path(s, path);
s = cd; s = [s, '/EFGChino']; path(s, path);
% clean memory and terminal
clear all
close all
clc

%% Boundary conditions
bc_homog = 1;  %non homogeneus (some nodes have BC imposed)
mockTheAgency = 0;

RHOx=3.2;

mesh = [7 11 21 41 81];
mesh_n = 4;
penalty = [10^3; (10^3)/8; (10^3)/64; 8*10^4; 10^4; (10^4)/8; 10^6; 10^7];
penalty_n = 8;
gaussOrderGL = 6;%cubature order: 2 (3 GPts), 4 (6 GPts), 6 (12 GPts)

condsK        = zeros(mesh_n, penalty_n);
ranges_n      = zeros(mesh_n, penalty_n);
errorsL2      = zeros(mesh_n, penalty_n);
numSolve      =  cell(mesh_n, penalty_n);
numGaussSolve =  cell(mesh_n, penalty_n);
analSolve     =  cell(mesh_n, penalty_n);

for imesh = 1 : mesh_n;
    for ipenalty = 1 : penalty_n;
        
    L       = [1 1];
    centre  = [0.5*L(1); 0.5*L(2)];
    Nx      = mesh(imesh);
    Ny      = mesh(imesh);
    h       = 1/(Nx-1);
    x_nodes = UniformGrid2D(Nx, Ny, L(1), L(2), centre);
    nPts    = size(x_nodes,1);
    
%% Gauss-Legendre Numerical Integration
% Numerical Integration
optGauss.orderGL     = gaussOrderGL;
bcOPTGaussY1.orderGL = gaussOrderGL;
bcOPTGaussY0.orderGL = gaussOrderGL;
bcOPTGaussX1.orderGL = gaussOrderGL;
bcOPTGaussX0.orderGL = gaussOrderGL;
% Sample points where the shape functions are computed
[x_gauss w_gauss] = MakeGLSamples2D(x_nodes, optGauss);

%% EFG paramters
optEFG.rho = RHOx*h;
optEFG.poliBaseType = '';
optEFG.order = 1;
% For compatibility------------------------------------------
% Options for the computation of the shape functions
optEFG.dim     = 2;           % spatial dimension
optEFG.verb    = 0;           % information and plots: 0=OFF, 1=ON
optEFG.grad    = 1;           % Computation of the Gradient 0:OFF 1:ON
optEFG.hess    = 0;           % Computation of the Hessian  0:OFF 1:ON
optEFG.TolNR   = 1.e-12;      % Newton-Raphson Tolerance
optEFG.Tol0    = 1.e-6;       % Target Zero Tolerance
% Definition of locality parameters
gamma  = 1.6;
beta   = gamma/h/h;
beta_n = ones(1,nPts)*beta;
optEFG.beta_n  = beta_n;      % locality parameters
range_n        = zeros(nPts,1);
for i=1:nPts
  range_n(i)  = max(2*h,sqrt(-log(optEFG.Tol0)./beta_n(i)));
end
optEFG.range_n = range_n;     % range definition
% end compatibility
%% for adjacency calculation, using four times the rho
range_n = ones(size(x_nodes,1))*4*optEFG.rho;

%%
% adjacency structure with the nearest neighbors nodes to each sample point
% nodal shape parameter
if (mockTheAgency == 0)
    s_near = SamplesAdjacency_efg(x_nodes, x_gauss, range_n);
else
    mockSnear = cell(1, size(x_gauss,1));
    for sx = 1 : size(x_gauss,1);
        mockSnear(sx) = {(1:size(x_nodes,1))};
    end;
    s_near=mockSnear;
end;
% Local-max entropy basis functions computation
optEFG.s_near = s_near;
outEFG = wrapperCHINO(x_nodes,x_gauss,optEFG);
p_samp  = outEFG.p_samp;
dp_samp = outEFG.dp_samp;

%% Modificacion para agregar puntos de gaus en el borde--------------------
% boundary Y=1-------------------------------------------------------------
count = 0;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,2) == 1)
        count = count + 1;
    end;
end;
bcNodesY1 = zeros(count,2);
i=1;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,2) == 1)
        bcNodesY1(i, : ) = x_nodes(ii, : );
        i = i+1;
    end;
end;
[bc_x_gauss1D bc_w_gaussY1]= MakeGLSamples1D(bcNodesY1(:,1), bcOPTGaussY1);
bc_x_gaussY1 = zeros(size(bc_x_gauss1D,1), 2);
bc_x_gaussY1(:,1) = bc_x_gauss1D(:);
bc_x_gaussY1(:,2) = 1;
% boundary Y=0-------------------------------------------------------------
count = 0;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,2) == 0)
        count = count + 1;
    end;
end;
bcNodesY0 = zeros(count,2);
i=1;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,2) == 0)
        bcNodesY0(i, : ) = x_nodes(ii, : );
        i = i+1;
    end;
end;
[bc_x_gauss1D bc_w_gaussY0] = MakeGLSamples1D(bcNodesY0(:,1), bcOPTGaussY0);
bc_x_gaussY0 = zeros(size(bc_x_gauss1D,1), 2);
bc_x_gaussY0(:,1) = bc_x_gauss1D(:);
bc_x_gaussY0(:,2) = 0;
% boundary X=0-------------------------------------------------------------
count = 0;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,1) == 0)
        count = count + 1;
    end;
end;
bcNodesX0 = zeros(count,2);
i=1;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,1) == 0)
        bcNodesX0(i, : ) = x_nodes(ii, : );
        i = i+1;
    end;
end;
[bc_x_gauss1D bc_w_gaussX0] = MakeGLSamples1D(bcNodesX0(:,2), bcOPTGaussX0);
bc_x_gaussX0 = zeros(size(bc_x_gauss1D,1), 2);
bc_x_gaussX0(:,2) = bc_x_gauss1D(:);
bc_x_gaussX0(:,1) = 0;

% boundary X=1-------------------------------------------------------------
count = 0;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,1) == 1)
        count = count + 1;
    end;
end;
bcNodesX1 = zeros(count,2);
i=1;
for ii = 1 : size(x_nodes,1);
    if (x_nodes(ii,1) == 1)
        bcNodesX1(i, : ) = x_nodes(ii, : );
        i = i+1;
    end;
end;
[bc_x_gauss1D bc_w_gaussX1] = MakeGLSamples1D(bcNodesX1(:,2), bcOPTGaussX1);
bc_x_gaussX1 = zeros(size(bc_x_gauss1D,1), 2);
bc_x_gaussX1(:,2) = bc_x_gauss1D(:);
bc_x_gaussX1(:,1) = 1;

nX0 = size(bc_x_gaussX0,1);
nX1 = size(bc_x_gaussX1,1);
nY0 = size(bc_x_gaussY0,1);
nY1 = size(bc_x_gaussY1,1);

nnX0 = size(bcNodesX0,1);
nnX1 = size(bcNodesX1,1);
nnY0 = size(bcNodesY0,1);
nnY1 = size(bcNodesY1,1);

% Put all the bcNodes all together
bc_x_gauss = zeros(nX0+nX1+nY0+nY1,2);
bc_w_gauss = zeros(size(bc_x_gauss,1),1);
bcNodes = zeros(nnX0+nnX1+nnY0+nnY1,2);

bcNodes(1:nnX0,:) = bcNodesX0(:,:);
bcNodes(nnX0+1:nnX0+nnX1,:) = bcNodesX1(:,:);
bcNodes(nnX0+nnX1+1:nnX0+nnX1+nnY0,:) = bcNodesY0(:,:);
bcNodes(nnX0+nnX1+nnY0+1:end,:) = bcNodesY1(:,:);

bc_x_gauss(1:nX0,:) = bc_x_gaussX0(:,:);
bc_x_gauss(nX0+1:nX0+nX1,:) = bc_x_gaussX1(:,:);
bc_x_gauss(nX0+nX1+1:nX0+nX1+nY0,:) = bc_x_gaussY0(:,:);
bc_x_gauss(nX0+nX1+nY0+1:end,:) = bc_x_gaussY1(:,:);

bc_w_gauss(1 : nX0) = bc_w_gaussX0(:);
bc_w_gauss(nX0+1 : nX0+nX1) = bc_w_gaussX1(:);
bc_w_gauss(nX0+nX1+1:nX0+nX1+nY0) = bc_w_gaussY0(:);
bc_w_gauss(nX0+nX1+nY0+1:end) = bc_w_gaussY1(:);

%% The boundary conditions and penalty values
bc = bc_x_gaussY1(:,1) - bc_x_gaussY1(:,1) .* bc_x_gaussY1(:,1);
boundaryConditions = zeros(nX0+nX1+nY0+nY1,1);
if (bc_homog == 1);
    boundaryConditions(nX0+nX1+nY0+1:end) = bc(:,1);
end;

%% Sampling the shape function on the boundary sample nodes----------------
% adjacency structure with the nearest neighbors nodes to each sample point
% nodal shape parameter
if (mockTheAgency == 0)
    bc_s_near = SamplesAdjacency_efg(x_nodes, bc_x_gauss, range_n);
else
    bc_mockSnear = cell(1, size(bc_x_gauss,1));
     for sx = 1 : size(bc_x_gauss,1);
         bc_mockSnear(sx) = {(1:size(x_nodes,1))};
     end;
     bc_s_near=bc_mockSnear;
end;
% Local-max entropy basis functions computation
bc_optEFG = optEFG;
bc_optEFG.s_near = bc_s_near;
bc_outEFG = wrapperCHINO(x_nodes, bc_x_gauss, bc_optEFG);
bc_p_samp  = bc_outEFG.p_samp;
bc_dp_samp = bc_outEFG.dp_samp;

%% The numerical solution is computed
%id_bd: indices for the nodes on the boundary
nPts    = length(x_nodes);
ids     = (1:nPts);
id_bd   = [ids(abs(x_nodes(:,1))<1.e-6) ids(abs(x_nodes(:,2))<1.e-6) ...
           ids(abs(x_nodes(:,1)-L(1))<1.e-6) ids(abs(x_nodes(:,2)-L(2))<1.e-6)];
id_bd   = unique(id_bd);

%% fix the h dependent penalty
        PENALTY=penalty(ipenalty);
        if (ipenalty==2)
            PENALTY=penalty(ipenalty)/h;
        end;
        if (ipenalty==3)
            PENALTY=penalty(ipenalty)/ (h^2);
        end;
        if (ipenalty==5)
            PENALTY=penalty(ipenalty)/h;
        end;
        if (ipenalty==6)
            PENALTY=penalty(ipenalty)/ (h^2);
        end;
% calculates the boundary nomals
% Caclulates the boundary normal at each sample point
boundaryNormals = zeros(nX0+nX1+nY0+nY1,2);
boundaryNormals(1:nX0,1) = -1;
boundaryNormals(1:nX0,2) = 0;
boundaryNormals(nX0+1:nX0+nX1,1) = 1;
boundaryNormals(nX0+1:nX0+nX1,2) = 0;
boundaryNormals(nX0+nX1+1:nX0+nX1+nY0,1) = 0;
boundaryNormals(nX0+nX1+1:nX0+nX1+nY0,2) = -1;
boundaryNormals(nX0+nX1+nY0+1:end,1) = 0;
boundaryNormals(nX0+nX1+nY0+1:end,2) = 1;

%% Standard solver with penalty boundary condition imposition
[u_nodes condK] = codinaSolveMatrix_BC_nitsche(bc_x_gauss, boundaryConditions, bc_p_samp, bc_w_gauss, bc_s_near, bc_dp_samp, boundaryNormals,...
                  PENALTY, bc_homog,nPts,id_bd,x_nodes,s_near,p_samp,dp_samp,x_gauss,w_gauss );

%% Error in L2 norm
% compute value of u in the Gauss points
u_gauss = SamplesSolution(u_nodes,s_near,p_samp);
% analytical value of u in the Gauss points
u_sol   = ex23_AnalyticalSolution(bc_homog,x_gauss);
% L_2 norm of the error
L2_GL   = sqrt(sum((u_gauss-u_sol).^2.*w_gauss));
fprintf('\tNorm L2 of the error : %10.4e\n', L2_GL);
figure(2)
plot3(x_nodes(:,1),x_nodes(:,2),u_nodes,'ro',...
    x_gauss(:,1),x_gauss(:,2),u_gauss, 'bx',...
    x_gauss(:,1),x_gauss(:,2),u_sol, 'k+')
xlabel('X')
ylabel('Y')
zlabel('u(x)')
legend('node points','Gauss points num','Gauss points real')

%% Solution is plot
figure(3);clf
nSY = 50;
nSX = 50;
plotLME2D(x_nodes,u_nodes, nSX, nSY, L(1), L(2), centre, optEFG);
view([20 10])
title(strcat('LME  \gamma = ',num2str(gamma)),'fontsize',18,'fontweight','b')
xlabel('X', 'fontsize',16,'fontweight','b')
ylabel('Y', 'fontsize',16,'fontweight','b')
zlabel('\delta_z','fontsize',16,'fontweight','b')
disp('- - - - - - - -')

        % salvar la solucion numerica y la analitica, salvar los errores
        errorsL2(imesh,ipenalty)       = L2_GL;
        numSolve{imesh,ipenalty}       = u_nodes;
        numGaussSolve{imesh,ipenalty}  = u_gauss;
        analSolve{imesh,ipenalty}      = u_sol;
        condsK(imesh,ipenalty)         = condK;
        ranges_n(imesh,ipenalty)       = range_n(1);
    end; % gamma iteration    
end; % mesh iteration


%% plot the errors curves.
figure;
axisH = mesh(1 : mesh_n);
axisH = axisH.^(-1);
errH1 = errorsL2(:,1);
errH2 = errorsL2(:,2);
errH3 = errorsL2(:,3);
loglog(axisH, errH1, 'or-'); hold on;
loglog(axisH, errH2, 'og-');
loglog(axisH, errH3, 'ob-');
loglog([0.01 0.1],[0.0001 0.01], 'k-');

%% save the important data
save('./FW_outputs/efg_nitsche.mat', 'errorsL2',...
        'numSolve',...
        'numGaussSolve',...
        'analSolve',...
        'condsK',...
        'ranges_n');